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ABSTRACT 

We present results of two simulations of the convection zone, obtained by solving the 
full hydrodynamic equations in a section of a spherical shell. The first simulation has 
cylindrical rotation contours (parallel to the rotation axis) and a strong meridional 
circulation, which traverses the entire depth. The second simulation has isorotation 
contours about mid-way between cylinders and cones, and a weak meridional circula- 
tion, concentrated in the uppermost part of the shell. 

We show that the solar differential rotation is directly related to a latitudinal 
entropy gradient, which pervades into the deep layers of the convection zone. We also 
offer an explanation of the angular velocity shear found at low latitudes near the 
top. A non-zero correlation between radial and zonal velocity fluctuations produces a 
significant Reynolds stress in that region. This constitutes a net transport of angular 
momentum inwards, which causes a slight modification of the overall structure of 
the differential rotation near the top. In essence, the thermodynamics controls the 
dynamics through the Taylor- Proudman momentum balance. The Reynolds stresses 
only become significant in the surface layers, where they generate a weak meridional 
circulation and an angular velocity 'bump'. 

Key words: Differential rotation, solar convection zone, compressible turbulence 



1 INTRODUCTION: OBSERVATIONS AND 
SIMULATIONS 

1.1 Observations 

In the outer 28 % by radius of the Sun, known as the Solar 
Convection Zone (SCZ), most of the vertical energy trans- 
port is by convection. Due to the combination of convection 
and rotation, the gaseous body exhibits differential rotation 
i.e. non-uniform angular velocity. Helioseismology observa- 
tions (Scherrer et al 1995, Libbrecht 1989) of the nearly 10 
million acoustic p-modes that leak from the interior into the 
atmosphere, have provided through frequency splittings, in- 
formation on the angular velocity distribution of the solar in- 
terior. The latest results suggest that the isorotation surfaces 
in the solar convection zone are cone-like (aligned radially), 
in disagreement with most numerical simulations, which 
tend to produce cyUndrical contours parallel to the rotation 
axis. Specifically, such observations have revealed that: the 
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angular velocity near the equator first increases then gently 
decreases with depth, at mid-latitudes it is almost constant 
with depth, while at high latitudes it increases with depth; 
on the surface of the Sun, the angular velocity increases from 
the poles (35 day period) to the equator (25 day period) and 
meridional circulation ve, is much weaker than differential 
rotation (zonal velocity ~ 2km/s,ii9 ~ 25m/s). 



1.2 Simulations 

1.2.1 Computational hurdles 

Modeling the SCZ is a formidable task. Observations sug- 
gest: that the SCZ is highly stratified, spanning about 19 
pressure scale heights in depth, with the Mach number 
(square of the ratio of flow velocity to sound speed) ap- 
proaching unity at the top; the Prandtl number (ratio of 
the timescales of thermal to viscous diffusion) is extremely 
small (10~^); the equation of state is complex; and finally, 
the motion is highly turbulent, revealed by the large number 
of scales ranging from 0.1 km to 10^ km (depth of SCZ). In 
the photosphere (near the top of the SCZ) the kinematic vis- 
cosity of hydrogen is about O.lcm^s"^, while the root mean 
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square (r.m.s.) velocity in the SCZ is the order of lOOms^^. 
These values suggest a Reynolds number, Re=velocity x 
length / kinematic viscosity, of at least 10^'^. As the num- 
ber of degrees of freedom needed to represent a flow is pro- 
portional to Re^''*, to resolve numerically the scales in all 
three directions would require 10^^ grid points. Even with 
present technology the maximum Reynolds number mod- 
eled in direct numerical simulations (DNS) is a few thou- 
sand (Brummel, Hurlbert & Toomre 1996). Alternatively, if 
one assumes most of the energy is contained in the resolved 
scales, which are much larger than the viscous and thermal 
dissipation scales (unresolved scales), and that the average 
energy transferred by these smaller scales can be modeled 
by entropy diffusion, then a much larger viscosity (smaller 
Re) can be used to model the turbulent flow. This is the 
principle of large eddy simulations (LES). 

1.2.2 Numerical simulations of the SCZ 

Explicit integration of the hydrodynamic equations requires 
the timestep to be less than the time for a signal (sound 
wave) to travel between two grid points. This is known the 
Courant-Friedrichs-Levy (CFL) stability criterion. Previous 
investigations have solved approximate forms of the fully 
compressible equations, in which sound waves do not exist. 
This circumvents the use of restrictively small timesteps. 
Incompressible convection in a shell was studied by Oilman 
(1978) using the Boussinesq approximation. This allows den- 
sity variation only when coupled with gravity. A better ap- 
proximation is to exclude only temporal variations in den- 
sity, known as the anelastic approximation. This still sup- 
presses sound waves but allows density to vary vertically. 
The most recent anelastic simulation in a full shell (global 
model) on massively parallel architectures are by Miesch et 
al (2000). Though there were signiflcant improvements on 
earlier work (Glatzmaier 1987), they were still unable to 
adequately resolve the observed cone-like structure of the 
angular velocity contours. Overall the isorotation contours 
were cylindrical. Furthermore, the angular velocity shear 
layer near the top of the convection zone could not be re- 
solved. Near the top, vigorous compressible motions have 
Mach numbers close to one and the anelastic approximation 
breaks down. 

The fully compressible equations have however been 
solved in local models (boxes). By considering only a tiny 
section of a spherical shell, Brummel et al. (1996) did a 
DNS of turbulent compressible convection in f-planes. In es- 
sentially the same geometry, Chan (1999) did a set of LES 
computations, in f-planes. The LES model emphasized ef- 
ficient convection (convective flux S> diffusive flux) as op- 
posed to the inefficient convection in the DNS models. In 
f-plane simulations the angle between the rotation vector 
and gravity is considered constant throughout the box. The 
problem with these local models is that they have periodic 
horizontal boundaries. In such domains the only possible 
source of differential rotation is the Reynolds stress. Hor- 
izontal averaging removes latitudinal gradients of thermal 
quantities which are a possible source of the differential ro- 
tation. For example, in the model by Durney (1999), the 
latitudinal entropy gradient (or baroclinicity) is essential in 
shaping the isorotation contours. Additionally, the f-plane 
models do not allow a realistic meridional circulation to de- 



velop; nor can Rossby waves exist as the rotation vector is 
constant throughout the domain. Rossby waves have been 
suggested as a possible source of correlation between longi- 
tudinal and meridional motions (Brummel et al. 1996). 

Our model represents a compromise between the anelas- 
tic global and compressible local simulations. We solve the 
full Navier Stokes equations in spherical coordinates in a 
domain with a significant latitudinal coverage. In this way 
we are able to resolve the upper shear layers and simul- 
taneously include the effects of meridional circulation and 
latitudinal gradients of thermal quantities. To circumvent 
the CFL restriction on the time step, semi-implicit time in- 
tegration is employed. We use the Alternating Direction Im- 
plicit Method on a Staggered Mesh (ADISM) (Chan & Wolff 
1982) in conjunction with an explicit method. Due to com- 
putational restraints of three-dimensional calculations, we 
are limited to studying a small section of the entire shell. A 
full simulation on a 70 x 70 x 39 three-dimensional mesh, 
spans just 60° in longitude and latitude. This calculation 
requires about a month to run when parallelized on 4 pro- 
cessors of the ORIGIN 2000. 

This paper consists of 4 more sections. Section 2 de- 
scribes the overall physical setup, formulates the mathemat- 
ical model and describes the numerical approach. This is 
followed by descriptions of the statistically averaged zonal 
and meridional flows, and turbulent nature of the compress- 
ible convection. The third section attempts to pin down the 
source of the differential rotation and meridional circulation. 
The flnal section is a conclusion. 



2 MATHEMATICAL MODEL 

2.1 Overall setup 

The computational domain is a 3-dimensional sector of a 
spherical shell, symmetric about the equatorial plane. A lon- 
gitudinal cross-section is shown in Fig. 0. The input energy 
flux at the bottom (straight arrows), is transported mostly 
by convective eddies (curly arrows) up to a height of 95% 
of the total depth. To pad out convective motions, a con- 
duction layer is placed in the upper 5% of the shell (shaded 
region). We consider two models that are identical apart 
from their longitudinal spans. The first model, denoted A, 
spans 30° in longitude and 60° in latitude, while the second, 
B, covers the same latitudinal range, but 60° in longitude. 
Both cover a total depth of 0.72R to R, which is the approx- 
imate depth occupied by the convection zone in the Sun. R 
is equivalent to the radius of the Sun. As depth is scaled by 
the total radius, the nondimensional radius is between 0.72 
and 1.0. Radial, meridional and zonal directions are labeled 
r (increasing outwards), 9 (increasing southwards) and (p 
(increasing eastwards), respectively. 

2.2 The equations 

In the absence of motion (d/dt = 0, v = 0), the equations 
governing conservation of mass, momentum and energy in a 
rotating stratified fiuid (Chan 1999), reduce to the equations 
of hydrostatic and thermal equilibrium. A solution of these, 
in spherical coordinates, is a polytrope 

T/Tt = l + Z{R/r~l)/{R/rb-l), (1) 
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Figure 1. Longitudinal cross-section of the computational shell 
with rotation axis f2. The input energy flux (straight arrows) is 
transported by convection (curly arrows) and then by radiation 
(shaded region), to the conducting top. 



p/p, = (T/TO", 
p/p, = {T/T,r+\ 



(2) 
(3) 



where r is the radial distance from the base, and n is the 
polytropic index. The subscripts 't' and 'b' denote a quan- 
tity measured at the top and bottom of the shell, respec- 
tively. T,p and p are the symbols for temperature, pressure 
and density. Z — (Tb — Tt)/Tt describes the extent of the 
stratification. 

Equations (|^^) provide a reference atmosphere from 
which appropriate dimensionless units can be formed. 
Length is sca led by the outer radius of the shell R , and 
time by R/ y^pt/pt- In such units, velocity is scaled by the 
isothermal sound speed at the top, y^pt/pt- From now on 
all quantities will be given in nondimensional units. Com- 
bining the equation of hydrostatic equilibrium, dp/dr = —pg 
and equation (jl|), gives in nondimensional units, gt = [n + 
l)Zrh/{l ~ Vb). We will consider gt as an independent pa- 
rameter. As the total number of pressure scale heights is 
(n -I- 1) ln(l -I- Z), the size of gt determines the depth of the 
layer. 

With such scaling, the governing equations become: 



dp/dt = 

dpv/dt — 

dE/dt = 

p = 



-V ■ pv, 

-V ■ pvv - Vp + V ■ S 
-pgr — 2pflo X V, 
1 



^ pTv + pv + (/Of /2)v 

-V • (-V ■ S + f) - pv ■ gr, 

pT, 



(4) 
(5) 



(6) 
(7) 



where f2o, S and f are defined below. Equations (^ 
represent a closed system of 5 dependent variables; density, 
radial meiss flux, meridional mass flux, zonal mass flux and 
total energy density, denoted by p, pVr, pve, pv^ and E, re- 
spectively. In a particular geometry (fixed r, 9 and (f) bound- 
aries), each turbulent convection simulation is specified by 
defining six nondimensional parameters. These are the ref- 
erence rotation rate Oo, the input energy flux at the base fh, 
the turbulent Prandtl number Pr = v/k, the gravitational 



acceleration at the top gt, the polytropic index n and the 
ratio of the specific heats 7. 

We will now consider individual terms in equations (^ 
Ignoring the coefficient of bulk viscosity (Becker 1968), 
the viscous stress tenser for a Newtonian fluid is Eij = 
IJb{dvi/dxj + dvj/dxi) — 2/i/3(V ■ 'v)Sij. In DNS, the viscos- 
ity is determined from the nondimensional parameters char- 
acterizing the convection simulation (e.g Pr, the Rayleigh 
number Ra, gtjTi and 7) . However, in LES, the dynamic 
viscosity p is increased so that it represents the effects of 
Reynolds stresses on the unresolved or sub-grid scales (SGS) 
(Smagorinsky 1963), 



p = p(c^A) (2fT : a- 



,1/2 



(8) 



The colon inside the brackets denotes tensor contraction of 



the rate of strain tensor ai 



+ VjVi)/2. 



The SGS 



eddy coefficient c^, is set to 0.2, the value for incompress- 
ible turbulence, and A'^ = rAr AS is an estimate of the local 
mesh size. Numerical experiments with different A and c^j 
are described in Chan & Wolff (1982). The present formu- 
lation ensures that the grid Reynolds number l\x v/v is oi 
order unity everywhere. To handle shocks, p is multiplied 
by 1 -I- 2/cJ'[{di.xdxVxY + (At/cJyVj,)^], where x and y de- 
note meridional and zonal directions (e.g. Aa; = rA6' and 
Ay = rsinOAij)) and Cs is the isothermal sound speed. As p 
is dependent on the horizontal divergence, any large hori- 
zontal velocity gradients are smoothed out by the increased 
viscosity. 

The gravitational acceleration at a distance r from the 
center of the sphere is gr, where g = gt/r^, and f is a 
unit vector directed radially outwards. The reference rota- 
tion vector is f2o = Q,od, where f2 is a unit vector in the 
direction of the rotation axis, i.e. Cl ■ r = cos9. The nondi- 
mensional rotation rate Qo, is equal to the ratio of the solar 
rotational velocity, to the isothermal sound speed at the top 
of the convection zone. For the sun Qo is about 0.3, so that 
the value of 2.91 used in the present computation, is associ- 
ated with rotational periods of about a day. 

The total energy per unit volume E, is the sum of the 
internal energy pT/(7 — 1) and the kinetic energy pv^ /2. 
The terms in the brackets on the left hand side of equation 
(^ are the various forms of energy flux in to or out of a 
unit volume fluid parcel. The first three terms constitute 
{E + p)v, which equals the convective fiux CppTw plus the 
kinetic energy flux pnP^ /2. Away from the upper and lower 
boundaries {E -\- p)v represents the energy transported by 
the resolved large scale eddies. The LES model is designed 
so that this term carries most of the vertical energy flux. 
The other fluxes from left to right, are the viscous flux and 
the diffusive flux, f . The last term in the energy equation is 
the work done by buoyancy. 

At the base f has a positive constant value. This is 
the source term of the vigorous turbulent convection i.e at 
r = Tb = 0.72, f = /br. At all other horizontal levels f acts 
as a diffusion term, computed as f = — fciVS — k2VT. The 
values of fci and k2 determine whether the layer is convective 
(unstable) or radiative (stable). In the unstable layer (0.72 < 
r < 0.986), fci = pT/Pr and k2 < /b/|VT|t=o. As k2 is very 
small, nearly all of the heat transport is convective. In the 
upper layer (0.986 <r< 1.0), fci = and fc2 = /b/|Vr|t=o, 
allowing conduction to transport all of the heat flux in the 
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stable layer. The conduction layer emulates radiation above 
the convection layer. 

The horizontal boundaries are insulating, stress free and 
impenetrable in latitude, and periodic in longitude. The 
top and bottom are both impenetrable and stress free. The 
source term /b is injected in at the bottom and conducted 
out through the top. This requires the top to be maintained 
at a constant temperature. 



2.3 LES versus DNS 

In the SCZ, as the convective flux S> diffusive flux, efficient 
mixing reduces the superadiabatic gradient V — Vad to just 
above zero (« 10~*). This type of convection is numerically 
simulated, by transferring energy from the base to the top of 
the computational domain as follows. Firstly, in the unstable 
layer (from r — 0.72 up to r = 0.986) the gas conductivity is 
artificially small, forcing convection to carry most of the heat 
flux across the imposed temperature gradient. Secondly, the 
initial polytrope is neutrally stable (7 = 5/3, n = 1.5). After 
the convection has developed, the relaxed thermal structure 
remains very close to this initial state, so that V is only 
slightly greater than Vad (or equivalently VS is just below 
zero). This approach to modeling the SCZ is very different 
from the DNS approach. The main difference lies in the role 
of the diffusive flux. In the DNS by Brummel et al. (1996) 
, radiation transports 75% of the energy flux, while the re- 
solved convective motions carry the rest (via the kinetic and 
enthalpy fluxes). In the present LES, the thermal conduc- 
tivity is nearly zero, therefore even the entropy diffusion of 
the SGS, is greater than radiative diffusion. As VS ~ 0, the 
resolved large eddies must carry the majority of the energy 
flux. 



2.4 Angular momentum conservation 

Due to the non-zero reference rotation rate, after relaxation 
one must ensure that there is no mean motion between the 
fluid in the shell and the rotating frame of reference. After 
the start of the computation, the compressible bulk can ac- 
quire some spurious form of rotation due to initial expansion 
or contraction. Enforcing the condition {pv^,)^ — 0, where 
'v' denotes volume averaging taken after thermal relaxation, 
ensures the total angular momentum is zero in the reference 
frame. This is accomplished by calculating the mean angu- 
lar velocity of the shell, (fi)^ = Epw^rsin^/Spr^sin^^, and 
subtracting the residual angular momentum from the total 
flow, pv^ — > pv,p — {Q)^pr^sm^9. When the statistics are 
gathered, (SI) ^ is less than 4 x 10"''. 



2.5 Numerical methods 

After some transformations to make the numerical scheme 
conservative and preserve second order accuracy for the 
nonuniform vertical grid (Chan & Sofla 1986), equations 
(^^) are discretised in spherical coordinates. Using a code 
developed by Chan and Wolff (1982), an implicit scheme 
(the Alternating Direction Implicit Method on a Staggered 
grid or ADISM) relaxes the fluid to a self consistent thermal 
equilibrium. The relaxation time is of the order of a 'Kelvin 
Helmholtz' time (« J pCyTdr/fh)- In the fully relaxed layer. 



the energy flux leaving the top of the shell is within 5 % of 
the input flux /b. 

Next a second order explicit method (Adams Bashforth 
time integration) gathers the statistics of the time averaged 
state. The statistical integration time is over 500 turn over 
times, and requires about 1.5 million time steps. For model A 
(longitudinal span of 30°), there are 39 x 70 x 35 grid points, 
in the radial , latitudinal and zonal directions respectively. 
For model B which has twice the longitudinal span of A, 
there are twice as many points in the longitudinal direction. 
The choice of grid comes from comparing a simulation in a 
very small section of the shell (Robinson 1999) {6 ± 15°), 
with a physically similar LES computation in a small box 
(Chan 1999). The box, which has an aspect ratio of 1.5, is 
placed at the mid-latitude of the shell. 

For a single processor on the ORIGIN2000, the CPU 
time per integration step is about 3 seconds in model B. Us- 
ing automatic parallelization on 4 processors, the speed up 
factor is about 3. Numerical stability requires a non-uniform 
grid with the same number of grid points per scale height 
and a very small input flux /b of 0.25/64. These signifi- 
cantly increase the total (implicit plus explicit) computation 
time. Consequently, the minimum time for a full simulation 
is about a month. 



3 RESULTS 
3.1 Trials 

After a series of numerical experiments (Robinson 1999) in 
a shell spanning ±15° in latitude and longitude, a set of 
parameters were found that generated a 'sun-like' rotation 
pattern. These were Qo = 2.91, /b = 0.25/64, Pr = 1 and 
gt = 19 (about 5 pressure scale heights). Fixing these param- 
eters, the latitudinal span was increased to 60° (30° above 
and below the equator). Surprisingly, the rotation profile 
appeared to be in some kind of 'quasi-steady state'. Ini- 
tially a 'sun-like' profile was seen, but as the computation 
progressed, the radial angular velocity gradient at the equa- 
tor, switched from positive to negative. This simulation was 
run to completion and the averaged results are classified as 
model A. 

Increasing Qo did not improve matters, and only when 
the longitudinal span was increased to 60°, was a positive 
angular velocity gradient sustained. The shell now covered 
the same longitudinal as latitudinal extent. This simula- 
tion was run independently of model A and the averaged 
results are classified as model B. We found out later that 
the flow reversal in A, was caused by a spurious meridional 
circulation. Strong downflows at the impenetrable latitudi- 
nal boundaries generate a powerful flow, pointed from the 
equ ator to wards the poles. This feature is described in sec- 
tion |3.2.2|. 



3.2 Flow characteristics 

In a turbulent fluid a quantity q can be split into a mean 
and a fluctuating part, 



q = q{r,e) + q'{r,e,<l,,t). 



(9) 
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The overbar represents a combined longitudinal and tempo- 
ral average, i.e. 

tl 

tl, is a time after the system has reached a self-consistent 
thermal equilibrium i.e. t\ > J edr/f^^. e is the internal en- 
ergy at each horizontal level. The time required for statisti- 
cal convergence {t2 — ti) depends on the particular quantity 
being averaged. While the mean velocities required about 
100 turnover times, turbulent quantities such as the veloc- 
ity correlation Vrv'^, took about 5 times as long. 

3.2.1 Mean zonal flow 

The angular velocity averaged over time and longitude, rel- 
ative to the rotating frame of reference, is coniputed and 
shown in Figs ^ and g for model A, and in Figs ^ and ^ for 
model B. In model A, the isorotation contours are parallel 
to the rotation axis. This means the angular velocity of a 
fluid element at any point in the shell is determined by its 
perpendicular distance from the rotation axis. Fig. |^ shows 
the angular velocity plotted against nondimensional depth 
(given as a fraction of the total solar radius, i.e. r = r/R) 
in the Northern hemisphere. Co-latitudes (90°-latitude) of 
90°, 85°, 79° and 67.5° are labelled by solid, long-dashed, 
triple-dot dashed and dot-dashed curves, respectively. The 
plots show that the angular velocity decreases radially out- 
wards and away from the equator, in direct contrast to the 
rotational profile found in the SCZ. 

If the simulation is repeated with twice the longitudi- 
nal span, the differential rotation has a much more 'sun-like' 
appearance. In model B, the shape of the isorotation con- 
tours resemble helioseismology observations in two distinc- 
tive ways. 

Firstly, an initial increase then decrease in angular ve- 
locity from the top inwards near the equator is found, as 
shown by the two closed circular rotation contours (or equiv- 
alently, the radial angular velocity profile in Fig. |^. He- 
lioseismic results, Kosovichev et al (1997), show that the 
increase inwards of angular velocity occurs at the equator 
and up to 60 degrees co-latitude, but the jury is still out 
on the behavior at higher co-latitudes, Schou et al (1999). 
As the computational shell only extends to about 60 de- 
grees co-latitude, we will only consider behavior away from 
the boundaries as being representative of the actual convec- 
tive flow. In that sense, within ±8° about the equator, the 
computed angular velocity 'bump' is at least qualitatively 
similar to the observed result. 

Secondly, away from the equator towards mid-latitudes, 
the contours are about half way between the cylindrical con- 
tours (Taylor Columns) seen in most earlier global simula- 
tions (e.g. Glatzmaier 1987), and the cone-like shape ob- 
served in the SCZ. There is also good agreement with recent 
anelastic global LES by Elhott, Miesch & Toomre (2000). 
Over the range of depth and latitude in common with their 
simulation and the present simulation, the isocontours are 
similar. The amount of variation of angular velocity with 
latitude is also in agreement with the SCZ. At the same 
colatitudes as in A, the mean angular velocity is plotted 



against depth. Fig. g. At the top of the shell the angular 
velocity, drops by about 0.2 between the equator (solid line) 
and co-latitude of 67.5° (dot-dashed line). As flo, is about 3, 
the drop implies a 7% variation in rotation rate over 22.5°, 
or extrapolating, a pole that spins 28% faster than the equa- 
tor. This is in rough agreement with the rotation rate at the 
surface of the Sun, which varies from 25 days at the equator 
to 35 days at the poles. 

3.2.2 Mean meridional flow 

The meridional flows are also very different. In A, a strong 
meridional circulation develops, directed from the equator 
towards the poles at the top, with strong downflows at the 
latitudinal boundaries. The mean meridional velocity vg at 
the same colatitudes as used for the zonal flow, is presented 
in Fig. ^ As W < in the Northern hemisphere, the upper 
leg of the circulation is directed from the equator to the pole. 
Moving polewards from the equator, where by symmetry the 
meridional velocity is zero, |W| increases to a maximum near 
a colatitude of 67.5°, and then reduces to zero as the flow 
approaches the latitudinal boundary. At the boundary, the 
magnitude of the maximum downward radial velocity (not 
shown) is about 0.2. The meridional circulation is produced 
by strong downflows associated with the stress free impene- 
trable boundaries. In early models of turbulent compressible 
convection in a small box, Chan & Sofia (1986), noticed that 
an impenetrable side boundary tends to attract downflows 
and they made all side boundaries periodic in later compu- 
tations. In a shell which excludes the poles, impenetrable 
boundaries are unavoidable. These strong downflows create 
an artiflcially large meridional circulation. 

Fortunately the downflows can be almost eliminated by 
widening the longitudinal span. While the magnitude of the 
maximum radial velocity at the boundary is 0.2 in A, it is 
less than 0.005 in B (both being directed downwards). The 
impenetrable boundaries have very little effect on the down- 
flows in B. This is because less restriction is placed on the 
flow direction. In A, the narrowness of the slice enhances 
the meridional circulation, enabling it to traverse the entire 
depth of the shell. In B, the meridional flow is only signifl- 
cant in the upper part of the shell (r > 0.95), elsewhere it is 
close to zero. Fig. |^ shows that ve is negative (poleward) at 
the top of the unstable layer, while just above, in the stable 
layer, there is a returning (equatorward) flow. The return 
flow is twice as fast as the poleward flow, because the fluid 
at the top has half the density of the equatorward moving 
fluid, i.e. momentum is conserved in the meridional velocity 
loop. The point where vg changes sign (r — 0.975) is very 
close to the beginning of the stable layer (r = 0.985). Be- 
low a depth of about about 0.95, W is very small and the 
multi-cellular appearance is most likely a residual effect. 

3.2.3 Turbulent quantities 

The root-mean-square (r.m.s.) variance of a quantity q is 
given by 

q" = (11) 

We can remove the radial dependence by averaging over the 
convection layer. 
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Figure 2. Isorotation contours averaged over time and longitude 
in a shell spanning 60° in latitude and 30° in longitude (model 
A). 



Figure 4. Isorotation contours averaged over time and longitude 
in a shell spanning 60° in latitude and longitude (model B). 



Figure 3. Depth variation of mean angular velocity in a shell 
spanning 60° in latitude and 30° in longitude (model A). Co- 
latitudes of 90°, 85°, 79° and 67.5° are denoted by solid, dashed, 
triple-dot-dashed and dot-dashed lines respectively. 



Figure 5. Depth variation of mean angular velocity in a shell 
spanning 60° in latitude and longitude (model B). Co-latitudes of 
90°, 85°, 79° and 67.5° are denoted by solid, dashed, triple-dot- 
dashed and dot-dashed lines respectively. 



(?"> = i y q"dr, (12) 

where d — 1 — is the total depth of the layer. The angled 
brackets will be used to denote depth averages. Note (g") 
depends only on latitude. 

Turbulent velocity characteristics at co-latitudes of 
90°, 84° and 67.5° are presented in tables |^ and |^ for 
the simulations A and B. The total kinetic energy (mean 

plus turbulent) is proportional to (v^) ^ , while the turbu- 
lent part of the kinetic energy is measured by {v"), where 

v" = u"^ + Vg^ + v'l^ . Columns 2 and 3 indicate the rel- 

1/2 

ative sizes of {v^) and {v"). In A, (u^) is more than 
4 times greater than (u"), i.e. most of the kinetic energy is 
in the mean circulation. Conversely in B, (w") is within a 

factor of 2 of {rP") , implying the kinetic energy is split 
about equally between the mean and turbulent scales. In 
model B, most of {v") is in the zonal component (f^'), while 
in A, {v'^) is the smallest component. The distribution of 
kinetic energy between the mean and turbulent flow may 
determine the nature of the differential rotation. By mea- 
suring the autocorrelation function of the vertical velocity, 
we found that the vertical scale of the turbulence is about 



1.5 pressure scale heights (P.S.H), compared to a total shell 
depth of 5 P.S.H. 

Nondimensional parameters defining the importance of 
rotation compared to other physical processes can be calcu- 
lated from the resultant flow. The strength of the turbulent 
convection relative to rotation is characterized by the Corio- 
lis number, Co — Qod/{v"}. The Reynolds number. Re, com- 
pares the relative magnitudes of the advection and viscous 
terms, and is calculated as iv")d/ Iji/'p). The importance of 
SGS viscosity relative to rotation, is measured by the Tay- 
lor number, Ta, which is AQ^^d"^ / (/l/p)^ For reference these 
non-dimensional parameters are computed from the relaxed 
flow and presented in the last three columns of the tables 
and ^. The Coriolis numbers are larger in A because the 
turbulent velocities are smaller. This suggests rotation will 
have a greater effect on the large scales in A than in B. The 
larger Reynolds numbers in B reflects the more turbulent 
nature of the flow. 



4 DISCUSSION: WHAT PRODUCES THE 
DIFFERENTIAL ROTATION? 
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Table 1. Model A: Turbulence characteristics in the 'non sun-like' differential rotation. 
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Figure 6. Depth variation of mean meridional velocity in a 
shell with a 30° longitudinal span (model A). Co-latitudes of 
90°, 85°, 79° and 67.5° are denoted by solid, dashed, triple-dot- 
dashed and dot-dashed lines respectively. 



Figure 7. Depth variation of mean meridional velocity in a 
shell with a 60° longitudinal span (model B). Co-latitudes of 
90°, 85°, 79° and 67.5° are denoted by solid, dashed, triple-dot- 
dashed and dot-dashed lines respectively. 



4.1 Taylor-Proudman balance (TPB) 

Averaging the meridional momentum equation over longi- 
tude and time, so that d/d<j) = and d/dt = 0, produces 



1 d , 2x 



d 



rsmO 
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Measured at an arbitrary co-latitude of 78°, Figs |8| and 
^ show the relative sizes of terms in (|l^) computed from 
the averaged flows in A and B. From left to right terms 
are labeled by pluses, stars, diamonds, triangles, boxes and 
crosses, respectively. The third term in each equation is split 
into two terms (denoted by the diamonds and triangles). For 
the Co, Re and Ta values in the present simulations, the vis- 
cous stress terms Ve are, the equator excluded, much smaller 
than other terms. Therefore these terms are not shown. 

In both A and B the latitudinal pressure gradient 
(boxes) and the Coriolis force (crosses) are in approximate 
balance. The terms which contain the Reynolds stresses are 
insignificant. The overall dominance of the pressure gradient 
and the Coriolis force suggest the Navier Stokes momentum 
equations can be approximated by: 



VP 



= g - X V. 



(14) 



The above equation is known as the Taylor-Proudman mo- 
mentum balance. Using entropy 5* = c„ln(p/p'') and the 
perfect gas equation p = pRT, the curl of the left hand side 
of (^^ produces: 



V X — = -\{VP X Vp) = VT X VS, 
p p^ 

so that the (j) component of the curl of equation 
alent to 



2f2o rcosO 



dr 



sinO 



do 



+ 



dr 89 86 dr 



(15) 



IS equiv- 



0. (16) 



As g can be written as a potential, the buoyancy term dis- 
appears. 

Measurements from the simulation show that dT/dr is 
nearly three orders of magnitude greater than 8T/d9. The 
temperature field is almost radial and does not alter much 
from the initial unperturbed state. Assuming T ~ T{r), the 
last term of ( |l6| ) can be neglected, and the second term can 
be replaced by 



dr 89 



9{r) dS 

Cn 89 



(17) 



which when substituted into (|16|), gives 
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V or oO 



(18) 



This is equivalent to the rotation law obtained by Durney 
(1999). 

Fig. hd shows the mean entropy, {S) averaged over time, 
longitude and depth, plotted against co-latitude for models 
A (triple dot dashed line) and B (solid line). In A, the en- 
tropy is almost independent of co-latitude, suggesting that 
dS/dd ~ 0, so that equation (|lq) reduces to 



dr 



sine^=0. 



(19) 



Equation ( |l9| ) implies v,f, = /(rsin6'), i.e. is constant along 
cylinders parallel to the rotation vector. In A because the 
pressure gradient and Coriolis terms dominate the momen- 
tum balance, and the latitudinal entropy gradient is approx- 
imately zero, Taylor columns are seen in the interior. This 
explains the isorotation contours in Fig. ^. 

Contrastingly, in B the entropy varies significantly with 
colatitude. This non-zero latitudinal entropy gradient is the 
reason why more 'sun-like' isorotation contours are seen in 
model B. In other words, the term —g{r)/cp ■ dS/d6 shapes 
the differential rotation. Following this line of reasoning, we 
computed the baroclinic term —g(r)/cp ■ dS/dO (denoted 
by boxes) and 2Qo{rcos6dv^ /dr — s\n6dv^/d9) (denoted by 
crosses) from the averaged flow. Fig. |l^ shows the results 
at co-latitudes of 90° (smallest magnitude), 82° and 72° 
(largest magnitude). Clearly they are in approximate bal- 
ance, confirming that the thermal structure i.e. Sir, 9), di- 
rectly controls the overall shape of the differential rotation 
profile. The Reynolds stresses play a only minor role in the 
dynamics. As the simulation is of 'mildly', rather than fully 
developed turbulence, it remains to be seen, whether the 
same is true in the sun. 



4.2 Meridional circulation and the Reynolds 

stress 

4.2.1 Meridional circulation produces differential rotation 
tn A 

In model A, meridional circulation drives the convection to 



a spurious equilibrium state (see sections 3T and 3.2.2 ) . The 
computation starts off with a sun-like rotational state, but as 
it progresses the layer relaxes to a completely different equi- 
librium. We previously suggested that this is a consequence 
of the artificially large downfiows at the impenetrable latitu- 
dinal boundaries. We will now show that the direction and 
strength of the meridional circulation in A is just enough to 
drive the differential rotation. 

Consider a fiuid parcel in a Lagrangian frame of ref- 
erence. Ignoring the zonal pressure gradient and frictional 
effects, the only force on the parcel is the Coriolis force 
associated with the meridional circulation {yr,vg,Q). The 
equation of motion of the fiuid parcel is then 



dv^/dt ~ ~2^1oCos9ve — 2^aSai9vr. 

Integrating w.r.t time gives 

Au^ ^ -2noA(sin6l) - 2noSin6IAr 



(20) 



(21) 



Figure 8. Successive terms in averaged meridional momentum 
equation measured at a co-latitude of 78° in a shell with a 30° 
longitudinal span (model A). Clearly the largest terms are the lat- 
itudinal pressure gradient (boxes) and the Coriolis force (crosses), 
while the Reynolds stresses are relatively insignificant. 



We can verify this relation by taking values of from Fig. 
h3 which shows how «^ varies with depth in model A. At the 



Figure 9. Successive terms in averaged meridional momentum 
equation measured at a co-latitude of 78° in a shell with a 60° lon- 
gitudinal span (model B). Similarly the largest terms are the lat- 
itudinal pressure gradient (boxes) and the Coriolis force (crosses) 
and again the Reynolds stresses are very small. 



top of shell (constant r), Au</,/A(sin&) ~ 5.9. While at the 
equator (constant 9), Au^/sinSAr ~ —5. As f2o = 2.91, the 
meridional circulation is about the right size and direction 
to produce the zonal velocity variation in A. This expresses 
as a tendency of fluid parcels to conserve their individual 
angular momentum. 



4.2.2 Reynolds stress produces meridional circulation in B 

In model B, neither the meridional circulation, nor the 
Reynolds stress, are sufficient to directly drive the differen- 
tial rotation. As previously described, it is the entropy dis- 
tribution which produces the differential rotation in B. The 
meridional circulation in B is completely different to that in 
A. Equation (^ij) cannot be applied to the radial and latitu- 
dinal variation of zonal velocity in B. Fluid parcels do not 
appear to conserve their individual angular momentum. In 
fact they do, it is just that other forces are involved, namely 
the Reynolds stresses. 

In a turbulent fluid, the velocity can be split into its 
mean and fluctuating part. Though the average of the fluc- 
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Figure 10. Longitude, time and depth averaged entropy varia- 
tion with co-latitude: simulations A and B are denoted by triple 
dot dashed and solid lines, respectively. The entropy is almost 
constant in model A, but has a significant latitudinal variation in 
B. In both cases the entropy gradient is zero at the equator. 



Figure 11. Averaged terms in the Taylor-Proudman balance for 
the shell with a 60° longitudinal span (model B). —g(r)/cp-dS/d9 
and 2Qo(''cos09D^/9r — sinddv^/ 89), are denoted by boxes and 
crosses respectively. Both expressions have the smallest magni- 
tude at the equator and successively larger values at co-latitudes 
of 82° and 72° . 



tuation is zero, the average of the product of two fluctuat- 
ing quantities is not necessarily zero. The Reynolds stress is 
the averaged correlation between small scale velocity fluc- 
tuations, in two component directions, and the density. To 
reduce the order of the statistical moments, we assume that 
density can be taken out of the correlation, 



(22) 

This approximation is good because the density fluctuations 
are small (of the order of the square of the turbulence Mach 
number) . This has been confirmed numerically. As the mean 
density only depends on depth, the nature of the Reynolds 
stress can be ascertained by looking at the velocity cor- 
relation, v'^v'j. Fig. |l^ shows v'rv'^ plotted at the same co- 
latitudes and using the same line markers as previous veloc- 
ity plots. Fig. |l| shows v'gv'^, at depths of 0.98, 0.95, 0.90 
and 0.85, denoted by solid, dotted, dashed and dot-dashed 
lines, respectively. 

By comparing the plots of (Fig. |^ and the veloc- 
ity correlations, certain features become apparent. Firstly, 
quantities are only significant near the top of the shell (above 
r = 0.90), elsewhere they are close to zero. Secondly, at lower 
co-latitudes (specifically 79° and 67.5°) the sign and order of 
magnitude of W closely matches negative the slope of vi-v'^. 
Thirdly, towards the equator, d / dr{v'rv'^) decreases, while 
d / d9{v'gv'^) becomes steeper (more negative). 

The relation between Ue and the Reynolds stresses (or 
velocity correlations) , can be traced to the zonal momentum 
balance. In a rotating frame of reference, the axisymmetric 
zonal momentum equation can be written as, 

9 ( ■ , 



Id. 2 
— — (i;r«</,r 
or 



N + 



rsin^ QQy"<'"'P — l + 
2Q,a{vgcos9 + Ursin^) ~ 0, 



(23) 



where A'' denotes the additional non-linear terms and the 
viscous terms are again excluded. Near the impenetrable 
top boundary, Vr is close to zero, so the last term can be ex- 
cluded. Clearly the derivatives of v'rv'^ and v'^v'^ are capable 
of driving Uq. Away from the equator v'rv'^ drives Ug, while 
close to the equator (specifically 85°) v'gv'^ becomes more 
important. 

Overall, the Reynolds stresses produce a weak merid- 
ional circulation concentrated in the uppermost part of the 
shell, but are negligible elsewhere. 



Figure 12. Longitude and time averaged zonal velocity in the 
narrower shell (model A). Co-latitudes of 90° , 85° , 79° and 67.5° 
are denoted by solid, dashed, triple-dot-dashed and dot-dashed 
lines respectively. 



4-2.3 Reynolds stress produces angular velocity 'hump' 
near the top 

One other feature of the zonal velocity profile that cannot be 
explained by the large scale interaction of convection with 
rotation (i.e. TPB), is the small angular velocity 'bump' 
(see Fig. 0). The initial increase in angular velocity moving 
inwards can be considered as a small scale feature of the flow 
(Au^ < 1% of the mean rotation rate). The 'bump' cannot 
be attributed to the large scale entropy variation. 

It seems likely that this small scale feature could be 
caused by the Reynolds stress. The plot of v'rv'^ (Fig- [l^ ), 
reveals a connection between a negative drop in v'rv'^ near 
the top of the sheU, and the zonal velocity variation (Fig. 

The negative trough near the equator at co-latitudes of 
90° and 85° (solid and dashed lines) represents an inward 
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Figure 13. Longitude and time average of the product of ra- 
dial and zonal velocity fluctuations v' rv' ^ (fs! Reynolds stress di- 
vided by mean density) versus depth (model B). Co-latitudes of 
90°, 85°, 79° and 67.5° are denoted by solid, dashed, triple-dot- 
dashed and dot-dashed lines respectively. The negative troughs 
near the equator (solid and dashed lines) represents an inward 
transport of angular momentum by the Reynolds stress. The pos- 
itive peaks away from the equator (triple-dot-dashed and dot- 
dashed lines) correspond to a change in the sign of the merid- 
ional velocity, indicating that IJg depends on the radial derivative 



Figure 14. Longitude and time average of the product of merid- 
ional and zonal velocity fluctuations v'gv'^ Reynolds stress 
divided by mean density) versus co-latitude (model B). Depths 
of 0.98, 0.95, 0.90 and 0.85 are denoted by solid, dotted, dashed 
and dot-dashed lines respectively. Near the top (solid and dot- 
ted lines) the latitudinal gradient of v'gv' ^ increases towards the 
equator. 

transport of angular momentum by the Reynolds stress. We 
suggest that this is the cause of the slight increase in angu- 
lar velocity (or 'bump') moving inwards from the top of the 
convection zone. Furthermore, the maximum angular veloc- 
ity at the equator, occurs at almost the same position as the 
minimum of VrV'.. 



5 CONCLUSION 

The aim of this work is to reproduce solar differential rota- 
tion, by solving the equations of hydrodynamics, in a section 
of a spherical shell. To achieve this, we have done two numer- 



ical simulations. The first simulation initially has an angular 
velocity that decreases inwards, in agreement with the solar 
case. However, as the computation progresses, the radial an- 
gular velocity gradient changes sign, and the rotation rate 
increases inwards. The mean rotational structure consists of 
cylindrical isorotation contours and a sub-rotating equator. 
The switch is due to an artificially large meridional circula- 
tion, which itself is a result of strong downflows that occur 
at the impenetrable side boundaries. 

In the second computation the longitudinal span is dou- 
bled, so that it now equals the latitudinal span. The wider 
shell has much weaker downflows and a more turbulent flow. 
Under these conditions the rotation profile remains in the 
'sun-like' state. The angular velocity now bears a closer re- 
semblance to solar case, with a radial and latitudinal vari- 
ation, both qualitatively and quantitatively similar to the 
SCZ. 

By using an implicit timesteping scheme, we are able 
to run the simulations for longer than a 'Kelvin Helmholtz' 
time scale. The emphasis on complete thermal relaxation, 
could be the vital ingredient, required to correctly model 
rotating convection. Incomplete thermal relaxation, may be 
one reason why earlier simulations, failed to get the proper 
rotation pattern. In the 'Sun-like' simulation, it is the large 
scale thermal structure (specifically the entropy variation 
with latitude), rather than small scale motions (as in mean 
field models or some numerical models), that directly pro- 
duces the differential rotation. The non-zero latitudinal en- 
tropy gradient (baroclinicity) shapes the differential rota- 
tion. We notice that in the 'sun-like' model, the kinetic 
energy distributes roughly equally between the mean and 
turbulent scales, in contrast to the non 'sun-like' model, in 
which nearly all of the energy is in the mean fiow. Further- 
more, the turbulent kinetic energy is greatest in the zonal 
direction. This suggests the turbulent nature of the SCZ, 
may have some indirect role in the maintenance of solar dif- 
ferential rotation, but how operates is still unclear. 

The Reynolds stresses are important in the upper 10 
% of the computational domain. This region contains the 
top of the convection (unstable) layer and all of the radia- 
tive (stable) layer. Here the velocity correlation vi-v'^, has a 
dual effect. Firstly, it provides an inward transport of an- 
gular momentum, which causes the slight increase in an- 
gular velocity, just below the surface and near the equa- 
tor. Secondly, —d / dr{v'rv'^) is principally responsible for the 
poleward meridional circulation at the top of the convec- 
tion layer, and the accompanying return flow in the stable 
layer. Closer to the equator (within about 5°), —d / dOiv'gv'^) 
appears to drive the meridional fiow. In the meridional mo- 
mentum equation, the equator excluded, the only significant 
terms are the pressure gradient and the Coriolois force. The 
Reynolds stresses are only important in the zonal momen- 
tum equation, here they generate the meridional circulation 
found at the top of the shell. 

The key result is that the dynamics near the top (sur- 
face layers) of the convection zone, are controlled by the 
Reynolds stresses, while elsewhere, the differential rotation 
is determined the Taylor-Proudman momentum balance. A 
similar conclusion was reached in the semi-analytical model 
of the SCZ by Durney (2000) (and references therein). 
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